module matrix
use strutils
use sysutils
implicit none 
!private 
!public  :: print_matrix 
!private :: print_matrix_real

!interface print_matrix
!   module procedure print_matrix_real
!end interface print_matrix

integer, private, parameter :: stdout_unit=6

contains

subroutine InitMatrix(xx, v)
   double precision, intent(inout) :: xx(:,:)   ! matrix to be printed
   double precision, intent(in) :: v
   integer ncol, nrow, i, j
   ncol = size(xx, dim = 1)
   nrow = size(xx, dim = 2)
   do i = 1, nrow 
      do j = 1, ncol 
         xx(j, i) = v
      end do
   end do
end subroutine InitMatrix

 
subroutine print_matrix(xx, out_unit, col_start, row_start, summarize, &
   col_width, num_dec)
implicit none 
! print a matrix of reals with optional integer and column labels and an optional title
   double precision             , intent(in)           :: xx(:,:)   ! matrix to be printed
   integer          , intent(in), optional :: out_unit      ! unit to which output printed --- standard output by default
   integer          , intent(in), optional :: col_start
   integer          , intent(in), optional :: row_start
   logical          , intent(in), optional :: summarize
   integer          , intent(in), optional :: col_width
   integer          , intent(in), optional :: num_dec
   character (len=100)                     :: fmt_row_use, fmt_col_use, fmt_total_use
   integer                                 :: i, j, nrow, ncol, out_unit_use, col_width_use, num_dec_use, &
                                                 ncolwithsum
   logical                                 :: summarize_use
   integer          :: col_start_internal
   integer          :: row_start_internal
   double precision :: f
   
   
   ncol = size(xx, dim = 1)
   nrow = size(xx, dim = 2)
   col_start_internal = 1
   row_start_internal = 1
   if (present(col_start)) col_start_internal = col_start
   if (present(row_start)) row_start_internal = row_start
   if (present(summarize)) then
      summarize_use = summarize
   else
      summarize_use = .false.
   endif
   if(present(col_width)) then
      col_width_use = col_width
   else
      col_width_use = 10
   endif
   if(present(num_dec)) then
      num_dec_use = num_dec
   else
      num_dec_use = 4
   endif
   
   ncolwithsum = ncol
   if(summarize_use) ncolwithsum = ncolwithsum + 1
   fmt_row_use = "(i6," // lrtrim(inttostr(ncolwithsum)) // "f" // lrtrim(inttostr(col_width_use)) // &
      "." // lrtrim(inttostr(num_dec_use)) // ")"
   
   fmt_col_use = "(6x," // lrtrim(inttostr(ncol)) // "i" // lrtrim(inttostr(col_width_use))
   if (summarize_use) then
      fmt_col_use = lrtrim(fmt_col_use) // ",A" // lrtrim(inttostr(col_width_use))
   endif
   fmt_col_use = lrtrim(fmt_col_use) // ")"

   fmt_total_use = "(A6," // lrtrim(inttostr(ncolwithsum)) // "f" // lrtrim(inttostr(col_width_use)) // &
      "." // lrtrim(inttostr(num_dec_use)) // ")"

   if (present(out_unit)) then
      out_unit_use = out_unit
   else
      out_unit_use = stdout_unit
   end if
   if (summarize_use) then
      write (out_unit_use, fmt_col_use) (i + col_start_internal - 1, i = 1, ncol), "SUM"
   else
      write (out_unit_use, fmt_col_use) (i + col_start_internal - 1, i = 1, ncol)
   endif
   
   do i = 1, nrow
      write(out_unit_use, '(i6)', advance='no') i + row_start_internal - 1
      do j = 1, ncol
         f = xx(j, i)
         call WriteFToFile(out_unit_use, f, col_width_use, num_dec_use, (j == ncol) .and. (.not. summarize_use))
      enddo
      if(summarize_use) then
         f = sum(xx(:, i))
         call WriteFToFile(out_unit_use, f, col_width_use, num_dec_use, .true.)
      endif
      
      !if (summarize_use) then
      !   write (out_unit_use,fmt_row_use) i + row_start_internal - 1, xx(:, i), sum(xx(:, i))
      !else
      !   write (out_unit_use,fmt_row_use) i + row_start_internal - 1, xx(:, i)
      !endif
   end do
   if (summarize_use) then
      write(out_unit_use, "(A6)", advance = 'no') "SUM"
      do j = 1, ncol + 1
         if (j == ncol + 1) then
            f = sum(xx(:,:))
         else
            f = sum(xx(j, :))
         endif
         call WriteFToFile(out_unit_use, f, col_width_use, num_dec_use, j == ncol + 1)
      enddo
   !   write (out_unit_use, fmt_total_use) "SUM", (sum(xx(i, :)), i = 1, ncol), sum(xx(:,:))
   endif
end subroutine print_matrix


end module matrix


!program xprint_matrix
!use print_matrix_mod, only: print_matrix
!implicit none
!integer          , parameter :: nrow = 5, ncol = 3, out_unit=10 
!character (len=*), parameter :: 
!fmt_title="(/,10x,'x')",fmt_col="(100i10)",fmt_row="(100f10.4)" 
!double precision                         :: xx(nrow,ncol) 
!logical          , parameter :: tf_row=.false. 
!call random_seed() 
!call random_number(xx) 
!! prints title, column and row numbers 
!call print_matrix(xx,fmt_title=fmt_title) 
!! prints title, row numbers 
!call print_matrix(xx,fmt_title=fmt_title,fmt_col="") 
!! print column numbers 
!call print_matrix(xx,fmt_col=fmt_col,fmt_row=fmt_row,print_row_int=tf_row) 
!! print title, row numbers 
!call print_matrix(xx,fmt_title=fmt_title,fmt_col=fmt_col, & 
!                  fmt_row=fmt_row,print_row_int=tf_row) 
!! print title, row numbers to unit out_unit 
!call print_matrix(xx,fmt_title=fmt_title,fmt_col=fmt_col, & 
!                  fmt_row=fmt_row,print_row_int=tf_row,out_unit=out_unit) 
!end program xprint_matrix 

